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I. INTRODUCTION 

With the experimental realization of Bose-Einstein 
condensation with neutral atomic gases |]-||], it has be- 
come possible to observe fundamental properties of quan- 
tum statistics directly. While there are numerous well- 
established fields of low-temperature quantum physics 
that deal with many-particle systems, most of the salient 
features related to the paradigm of indistinguishablility 
are masked by strong interactions, or contact with ther- 
malizing reservoirs. The newly gained ability to isolate 
the condensate fraction from its environment to a high 
degree, as well as the tight spatial confinement of the 
macroscopic quantum field, has opened up ways to selec- 
tively manipulate the mean-field, and study its temporal 
evolution towards equilibrium almost " in vivo" . 

The motion of trapped atoms in a dilute gas consists 
of free oscillations within the external potential that are 
interrupted by short binary collision events. Convention- 
ally, that interaction strength is measured by the range 
of the repulsive two-particle potential, i.e. the scattering 
length a s . The duration of a collision event To is given by 
the time a particle of average velocity v spends in the in- 
teraction region, i.e., To = a s /v. On the other hand, the 
inverse collision rate r c = l/(na^ v) is an estimate of the 
time between successive collisions where n denotes a par- 
ticle density. As we are interested in low kinetic energies 
and the weak interaction limit, one finds a characteristic 
separation of time scales, i.e., 



r c > r , or na s <C 1. 



(1) 



It is this separation of time-scales that gives raise to a 
kinetic stage of evolution, preceding any equilibrium sit- 
uation. For example, such a kinetic stage is absent in 
strongly interacting systems where a local equilibrium is 
established immediately (r c w To, hydro-dynamic stage). 
Any individual collision event creates a quantum mechan- 
ical entanglement between collision partners. However, 
due to the long separation between successive collisions 
and the presence of intermediate weak fluctuations these 



temporal correlations decay rapidly (Markov approxima- 
tion). 

Based on this microscopic picture of the weakly inter- 
acting bosonic gas, we derive a generalized kinetic theory 
for a coarse-grained Markovian many-particle density op- 
erator, as discussed by A. I. Akhiezer and S. V. Pelet- 
minskii, [Q. The quint-essential assumption behind this 
approximation is a coarse-grained density operator that 
depends only on a few selected variables. These quan- 
tities will serve as master variables and determine the 
system's evolution on a coarse-grained time scale. From 
a perturbative expansion of this density operator, we ob- 
tain kinetic equations that describe the temporal evolu- 
tion of the expectation value of any (single-time) observ- 
able. As a specific application of this formulation, we 
then assume that the condensed gas can be described es- 
sentially by the dynamic evolution of mean fields, normal- 
and anomalous fluctuations. In particular, we obtain ki- 
netic equations that comprise the Gross-Pitaevskii equa- 
tion and the quantum-Boltzmann equation as special in- 
stances. These kinetic equations are formulated basis- 
independently and consider all second-order processes 
which give raise to collisional energy shifts and damp- 
ing rates. It is also noteworthy that the present theory 
can be applied readily to multi-component bosonic gases, 
provided the two-particle scattering matrix is interpreted 
accordingly. 

However, considering the rapid development of this 
field of physics, it is also of great importance to compare 
the predictions of this and other kinetic theories |^ [p]] 
(mean-field equations and kinetic theories) with the re- 
sults obtained from finite-temperature calculations based 
on the Hartree-Fock-Bogoliubov formulation |jp6|- p5| (col- 
lective excitation frequencies and damping rates), direct 
approaches to solve the many-particle Schrodinger equa- 
tion |26|27j] (configuration interaction, hyper-spherical 
coordinates), renormalization group techniques |B8fl, and, 
above all, to gauge them against physical reality |p9|-pl[ . 
Further references on the extensive literature can be 
found in recent review articles and literature compila- 
tions isyji. 



The present article is organized as follows: In Sec. [I A 



we introduce a coarse-grained statistical density opera- 
tor and derive an integral equation for it in Sec. II B 



Assuming a weak two-particle collision rate, permits the 
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develop ment of a perturbation theory, which is described 
in Sec. II C . From an explicit series expansion of the 



II D 



coarse-grained density operator, we obtain in Sec 
a set of kinetic equat ions c haracteristic of a condensed 
bosonic gas. In Sec. Ill A, we introduce the Hamilton 



operator that governs the kinetic evolution of a weakly 
interacting, repulsiv e gas . A set of relevant operators 



is introduced in Sec. Ill B , i.e., macroscopic mean-fields, 



normal fluctuations, as well as anomalous averages. By 
applying the general kinetic master equations to this spe- 
cific set of relevant operators, we obtain self-consistent ki- 
netic master-equations and give a detailed discussion in 
Sec. IIIC. Finally, in Sec. [V, we discuss the prospects 



and possible implications for a numerical simulation of a 
fully quantum mechanical self-consistent solution. 



II. A COARSE-GRAINED MARKOVIAN 
DENSITY OPERATOR 

The quantum mechanical state of a many-particle sys- 
tem contains an overwhelming amount of data and it 
is the objective of quantum statistical mechanics to ex- 
tract the relevant information from a data subset as small 
as possible. Such a minimal description then leads to a 
coarse-grained picture of the behavior of the system. The 
effects of coarse-graining on the dynamical evolution are 
well known from the quantum mechanics of open sys- 
tems ]3^|36|]. By discarding some observable informa- 
tion from the unitary evolution of the complete system, 
one breaks the time-reversal symmetry, thus introducing 
irreversibility. However, a judicious partitioning of the 
many, interacting degrees of freedom of a large system 
into a small relevant subsystem and a weakly coupled, 
complementary reservoir leads to a tractable description 
of the subsystem's evolution towards equilibrium. This 
equilibrium state is determined by the properties of the 
reservoir. 

In this section, we will pursue these general ideas and 
define a coarse-grained statistical density operator that 
depends functionally on only a few fundamental vari- 
ables, and thus gives a raw picture of the "true" state 
of the many-particle system. Based on the definition for 
the coarse-grained statistical density operator, we then 
derive an integral equation that determines the functional 
form of this operator. In the limit of weakly interacting, 
dilute quantum gases where strong collisional interaction 
events are well separated in time, it is possible to solve 
this equation perturbatively and establish a hierarchy in 
terms of an expansion parameter proportional to this in- 
teraction strength. Once a perturbative expression for 
the functional form of the coarse-grained statistical op- 
erator is determined, we then use it to study the motion 
of expectation values of general observables. These equa- 
tions of motion establish a closed, self-consistent set of 
kinetic equations for the following restricted set of fun- 
damental variables: the mean field, the normal single- 
particle density, and the anomalous fluctuations. By this 



method, we generalize the Gross-Pitaevskii mean-field 
equation and merge it consistently with an approach that 
leads to the quantum-Boltzmann equations for both nor- 
mal densities and anomalous fluctuations. 

In many ways our approach is reminiscent of the clas- 
sical Bogoliubov-Born-Green-Kirkwood-Yvon (BBYGK) 
method | ]37fl , for in a similar manner, higher order correla- 
tions can also be characterized by a restricted set of vari- 
ables. Specifically, in classical Markovian systems higher 
order correlation functions can be expressed in terms of 
single-particle densities. 



A. Fundamental assumptions of statistical mechanics 

The derivation of a coarse-grained density operator re- 
lies on the basic assumption of statistical mechanics that 
any non-equilibrium statistical correlation eventually de- 
cays P,p^7|-^0|. For example, in classical kinetics, this 
time is typically of the order of the duration of a colli- 
sion. 

In the following discussion, we introduce three distinct 
many particle density operators. The first density op- 
erator, p(t), describes the complete state of the many- 
particle system. Starting from an initial value p(0), the 
system evolves unitarily in time according to a time in- 
dependent Hamilton operator H. In analogy with the 
celebrated Chapman-Enskog |3^,|3£| procedure of classi- 
cal statistical mechanics which was later introduced to 
quantum statistics by Bogoliubov |Q] , we assume that all 
characteristic features of the quantum statistical state 
p(t) can be inferred from expectation values of a certain, 
restricted set of operators {7i|? S I}, where the index 
i enumerates linearly independent operators ji from a 
(possibly infinite) index set X. 

Thus, we introduce a second coarse-grained, Marko- 
vian density operator <T{ 7i (t)| iGi} that approximates the 
complete density operator p(t) 



-iHt 



p(0)e 



HI 



(2) 



The set of expectations values {7i(i)| i £ 1} that parame- 
terize the coarse-grained density operator are found by a 
quantum average (...) over all states of the many-particle 
configuration space 



(3) 



7t(t) = (7i) = Tr {% P^)} = Tr {% °"{7i(t)l jei}} • 

While the coarse-graining assumption, i.e., the restric- 
tion to a few selected variables, is a statistical statement, 
the Markovian postulate concerns the separation of the 
time scales that govern the processes of equilibration and 
the decorrelation of fluctuations. Within this limit, the 
temporal evolution of this coarse-grained statistical den- 
sity operator o^mitez} i s solely governed by the mo- 
tion of the set of expectation values {7i(i)| i £ 2T}. Any 
time dependence of the matrix elements of C{ 7 (*)} can be 
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attributed to the evolution of this restricted set of oper- 
ators {7^}, thus no relevant intrinsic time dependence is 
unaccounted for R. Although, the set of operators is un- 
specified so far, it can include operators such as unity (1), 
number (N), linear momentum (P), angular momentum 
(L), and energy (H). The larger the set of operators, the 
better the quality of any subsequent approximation. 

The third relevant many-particle density operator is 
(Tr°\. It serves as a reference distribution and describes 
a relaxed (but non-equilibrium) state of the gas between 
consecutive collision events. We define it to yield the 
same expectation values as cr{ 7 } on the restricted set of 
operators {74}. Thus, this reference distribution is given 
by 

*fj } = exp^r^), (4) 

where repeated indices imply a summation over operators 
7, and their conjugate thermo-dynamic coordinates T \ 
These coordinates are defined implicitly by 

7i = Tr {li ^{7} } = Tr {"Ji } ■ (5) 

It is of interest to note that the reference state, as given 
by the exponentiated form in Eq. (Q), is still of the most 
general form permitted by the principles of quantum me- 
chanics, as long as the set of operators {-ji} is complete. 

This ansatz for a non-equilibrium coarse-grained sta- 
tistical operator <j{ 7 } and a "self-adjusting" reference 

distribution 07 \ is also used in various contexts of 
{7} 

equilibrium thermo-dynamics. For example, from a fi- 
nite order truncation of a quantum virial expansion 
j|, one obtains a coarse-grained approximation cr^^py 
of the grand-canonical statistical density operator p = 
exp (fi — fiN — f3H). The starting point of the itera- 
tion procedure is a reference distribution tr|J| ^ that 
matches the solution as closely as possible. In here, the 
conjugate thermo-dynamic coordinates {f2, /z, /?} corre- 
spond to the observables: unity, number and total en- 
ergy. 

B. Derivation of an integral equation 

To find the functional form of the postulated coarse- 
grained Markovian density operator &{<y(t)}i we use a 
"boot-strapping" method. First, by assuming a given, 
yet unknown solution of Liouvillc's equation 



*In the absence of ambiguities, we simplify the notation of 
the set of operators or their corresponding expectation values 
by dropping the index set X, the indexing label i, or even the 
time argument t, completely. 



■J t a {i(t)} = -i[H, cr {7 ( t )}], (6) 

we can determine the dynamical evolution of the expec- 
tation values 7i(t) = Tr{7^ C{ 7 ( t )}} from 

j t7i (t)=iTr{[H,%]a Mt)} }. (7) 

Second, to derive the functional form of cr{ 7 (t)}j we now 
use Eq. (||), Eq. (0), and the fact that there is no explicit 
time dependence in C{ 7 (t)}- Thus, by partial differentia- 
tion one obtains the following equation 

Tr I [H, 7l ] cr {7} I d 1% a {l} = - [H , ] , (8) 

where again, we have adopted the convention that re- 
peated indices imply summation, unless stated otherwise. 
Moreover, by dropping the explicit time dependence of 
the reference point in phase-space {ji(t)}, one can con- 
sider this as a partial-differential equation with indepen- 
dent variables {"fi\- 

The total Hamilton operator H that governs the evo- 
lution of a weakly interacting, dilute gas permits a par- 
titioning of the energy into a free part and a pre- 
sumably weak interaction if' 1 ' 

H = £T (0) + H {1) . (9) 

One could use the bare, single-particle energy H^> that 
determines the free kinetic evolution of the gas as a start- 
ing point of a series expansion of the coarse-grained den- 
sity operator in terms of the interaction strength. How- 
ever, it is well known that the mean-field interaction will 
significantly affect the single-particle energies. Anticipat- 
ing this, we will shift the expansion point by an as 
yet undetermined single-particle energy Q^jy to a dressed 
energy 

(10) 

To conserve energy, we have to reduce the interaction 
energy by an equal amount 

HV=HV- Q«. (11) 

The shifted interaction energy represents the fluc- 

tuations about the mean-field energy and will be consid- 
ered as "weak", or in other words, that strong fluctua- 
tions are well separated in time. This procedure is anal- 
ogous to the first order energy shift found with ordinary 
Schrodinger-Rayleigh perturbation theory. The explicit 
form of this single-particle renormalization energy Q9^-i 
which is of the same order as the interaction energy, will 
be determined in the course of this calculation. 

By expanding Eq. (||) around the dressed single- 
particle energy Hi°\, one obtains 
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Tr 



%}a { , } } 9 7 , CT{7} + [H^ } , a h} ] = ^J, (12) 

ere Fr , is introduced for conv 

{7} 

maining first order contributions 



where F^l is introduced for convenience to hold the re- 

{7} 



final value, it also evolves freely according to Eq. (|lE|). 
To account for this evolution, we also define a unitary 
propagator ui°l (r) by 



^ ^ <t W} ] -Tr{[H$ } , %]<t {j} } d 7l a {l} . (13) 



{7} 



{7F 

d m°\(r) = -im , m fjf\{r), 

{ 7 }V J {7(t;{7})} h 



There are two interesting points in considering Eq. (12). 
First, there is not one unique solution cr{ 7 }, but solution 
manifolds that are labeled by the constants of motion. 
Thus, any particular solution has to be augmented with 
appropriate physical boundary conditions. Second, the 
structure of the commutator of the single-particle Hamil- 
tonian with the set of operators {%} is an intrinsic 
property of the system. It is determined both by the 
particular physical configuration and the set of opera- 
tors. If the set of operators is chosen appropriately, the 
commutator forms an algebra with structure constants 
^4{ 7 } defined by 



(14) 



If we substitute this algebraic closure relation into 
Eq. (O), one obtains 



(15) 



By the method of characteristics, one can transform 
this first order, inhomogeneous partial-differential equa- 
tion into an equivalent integral equation. Although for- 
mally equivalent, the later method is advantageous as 
it leads naturally to a series expansion of the coarse- 
grained density operator by iteration (compare Dyson se- 
ries pi]]). The characteristic trajectories are parameter- 
ized curves { 7i (r; {7})} along which the boundary value 
of the coarse-grained density operator is propagating in 
phase-space according to Eq. (|l^). They are defined by 

d 



X7i(T;{7» = i " 4 {7(r;{7})}/Ti(^{7}) I (16) 

(17) 



dr 

7,-(t = 0;{7}) = 7 



and the corresponding boundary conditions. Formally, 
the solution of the differential equation defines a map 
Kiy\(r) from regions in phase-space connected by the 
characteristics, i.e., 



(18) 



This map is obtained from the equation defining the char- 
acteristics Eq. ([L6J), i.e., 



(19) 



— K {l} (r) t J = iA mT , {7})} . 1 K {i} (t)i j 
K m (t = 0) = 1. 



(20) 



As the coarse-grained density operator is transported 
along the characteristics from its boundary value to its 



dr 
>(o) 



tM(t = 0) = 1. 



(21) 
(22) 



The solution for the inhomogeneous partial-differential 
equation Eq. ( |l5| ) is now obtained from an interaction- 
picture representation of the density operator u^{t) 
that is defined by 



^{7}M = U$(t) <x {7(t;{7})} U™(t). (23) 

At t = 0, it coincides with the Schrodinger-picture value 

a {7} (r = 0)=a {7} , (24) 

and, if our reference distribution, as defined in Eq. (^) 
matches the input state in the remote past asymptoti- 
cally, i.e., 

,(°) 



r(o) 



i!™ °{7(t;{7})} 



lim a,-, , 1 , 1 . 

^_oo {7(t;{7})}' 



then we find for the initial condition that 

,(0) 



Mm o- {7} (r) = o-y } . 



(25) 



(26) 



Finally, by differentiating the interaction-picture repre- 
sentation Eq. ( p3| ) , the use of Eq. ([l5|) and a subsequent 
integration subject to the boundary conditions given in 
Eqs. (^4 26), the integral equation for the coarse-grained 
density operator is obtained: 



(0) 



<7{ 7 } = (Tr.,! — I 



{7} 



I 

l7(r;{ 7 }) 



+ Tr{[^ } , 74 ]a {7} } d^^U^r). 



(27) 



This expression links the interacting coarse-grained den- 
sity operator to its non-interacting value by evaluating 
the integrand along its collision history. A regularizing 
function (77 — > 0+) has been introduced and suppresses 
higher order correlations, built up in individual collisional 
events. 



C. The limit of weakly interacting, dilute gases 

In the case of a weakly interacting system, we can seek 
the solution to the integral equation Eq. (|27j ) in form 
a power series of the density operator. The expansion 
parameter is given by the interaction strength 



ff {7> =E CT {tV 
1=0 



(28) 
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If this series is inserted into the integral equation for the 
coarse-grained density operator, then one finds for the 
first order correction 



r (1) 



i I dr e 

' — oo 



Tr 



" T ff$V) (29) 

f rS-W ~ i „(o) \ * jo) \ ff(o) , s 



We will show that the first two terms of the series Eq. (|2£ 
are sufficient to determine the kinetic equations for the 
expectation values of the set of operators {71} including 
second order (collisional) contributions. 

It is important to note that the integrand has to be 
evaluated along the trajectories to the past. However, it 
can be further simplified by using again the interaction- 
picture representation, as derived from Eqs. ( |l4| , |i9| , pl"| ) . 
First, one notices that the interaction-picture representa- 
tion of the reference distribution is constant, if evaluated 
along the characteristic trajectories. This is implied in 
its definition. 



AO) _ f)(o)t/ s (0) u {0) (t) 



(30) 



Second, one finds that in the interaction-picture the rel- 
evant operator ji(T;{'y}) can be expressed as a linear 
combination of Schrodinger-picture operators 7,- 

Ur; {7}) = U${t) 7, ^(r) = K {l} (r)^ % (31) 
Third, we define an interaction-picture Hamiltonian by 



{7} i v {7(r;{7 })}^ {7 ^ 



(32) 



With these definitions, one can evaluate the integrand 
along the characteristic. Details of this calculation are 
outlined in Appendix |A|. Thus, within the limit of weak 
interactions and the Markov approximation, one obtains 
the first two contributions to the coarse-grained density 
operator as 



(0) 



+ Tr{[H$ } (r), 7i]<\} d^a^) +0[2]. (33) 

It is easy to show that within these approximations the 
coarse-grained density operator is Hermitian, and that it 
yields exactly the same expectation values as the refer- 
ence distribution alone, as defined in Eq. (|^). 



D. Quantum kinetic equations 

With this first order result for the coarse-grained den- 
sity operator, we can now address the second part of 
our "boot-strapping" procedure: the kinetic evolution of 
an expectation value (6) subject to this coarse-grained 



density operator. The kinetic equations are obtained by 
averaging Heisenberg's equation with the coarse-grained 
density operator, as given in Eq. (|t]). As we are inter- 
ested in a power series expansion of the kinetic equations, 
again we decompose the total Hamilton operator H and 
the coarse-grained density operator o-{ 7 } into its various 
contributions in terms of the interaction strength, i.e., 



dt 



(6) = iTr{[H, o]cr {7(t)} | 



(34) 



By grouping the individual terms, one finds 



(35) 



1=0 



where we have introduced linear Liouville operators 



(36) 



L 



^) [6]=i Tr{[H^ } ,6]^}. (37) 



1. General observables 

To obtain practically applicable approximations for the 
kinetic equation, we will truncate these series at low or- 
der, i.e., first or second order. In the case of a general 
operator 6 that can not be represented by a linear com- 
bination of relevant operators, i.e., 6 (f. Span({7i | i G I}) 
this means that 



(o) 

Mm 



[d]+L 



(1) r 



(38) 



R M t)}[o]+L 



(2) r 



0[2], 



which is correct up to first order. At first glance, it seems 
to be inconsistent to include also the second order con- 
tribution L^ 2 ) in this first order expression. However, a 
closer inspection shows that this particular kinetic equa- 
tion Eq. ([38]) preserves all constants of motion. In other 
words, if the operator 6 is an exact symmetry 



then, according to Eq 
served 



[H, 0} = 0, (39) 
1), all initial averages are con- 



5< 6 > = a 



(40) 



This is particularly interesting if we consider constants 
of motion related to number (N,N 2 , . . .) and energy 
(H, H 2 , . . . ), as implemented in many recent BEC exper- 
iments. If the systems are prepared initially in one of the 
standard thermo-dynamical ensembles: micro-canonical 
(AN = 0,AE = 0), canonical (AN = 0, AE), or grand- 
canonical (AN, AE), then these properties are preserved 
in time. 
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2. Master variables 



In the previous section the kinetic evolution of a gen- 
eral operator was examined. Now, we focus on the set 
of relevant operators {7i}. The kinetic equations for the 
corresponding expectation values {7^} constitute a self- 
consistent set of master-equations that determine the sys- 
tem's evolution. Again, we are only interested in a low 
order truncation of Eq. (^5|). Due to the requirement 
that the reference distribution yields the same expec- 
tation values as the unexpanded state Eqs. (^,0), all 
contributions of [li] = vanish identically. Thus, 

we find a quantum kinetic master-equation correct up to 
third order, i.e., 

= *$t»fo] + L {%)}^ + L Mt»M + 0[3]. 

(41) 

By using the explicit expression for the coarse-grained 
density matrix Eq. (|33|), we find the following quantum- 
Boltzmann equation 

jp(t)=iTv{[H,^]a^ m } + (42) 
; dre^ Tr {<#> (t)} [^ (t)} (r), { %(O), + % (^L^jft] + Tr {[^#£^(0), % } <$ (t)} })] } , 

where _ffj*j(r) is the interaction picture Hamiltonian 
as defined in Eq. (|32]). 

There are two interesting features in Eq. ([42|): First, 
the coherent part depends only on the total Hamiltonian 
H. Eventually, this part will determine the evolution of 
the expectation values subject to the external trapping-, 
and the mean-field potentials. Consequently, all derived 
mean-field potentials are invariant with respect to a par- 
ticular partitioning of the total energy H , as introduced 
by the renormalization potential Qsj\ (see. Eq. (|l0|)). 
Second, the kinetic equation is strictly local in time, i.e., 
Markovian. The validity criterion for this result is a cor- 
relation time of the energy fluctuations which is much 
shorter than the time scale upon which the expectation 
values evolve and this is the case for a weakly interacting 
dilute gas. 
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III. KINETIC EQUATIONS FOR MEAN FIELDS, 
NORMAL AND ANOMALOUS FLUCTUATIONS 

In the previous sections, we determined the form of a 
coarse-grained density operator as a functional of a cer- 
tain set of mean values. In turn, the temporal evolution 
of the mean values is determined by the kinetic equations 
Eq. ([42l ) . In this section, we apply these general results to 
the particular situation of weakly interacting, low tem- 
perature bosonic gas. Thus, we have to determine the 
predominant energy contributions to the system Hamil- 
tonian H , as well as, decide on a relevant set of operators 

A. The dynamical evolution 

In second quantization, the removal or addition of a 
particle from or to a position x is described by the action 
of quantum field operators d x , a£. on the corresponding 
quantum state. As we are considering bosonic particles, 
these fields have to obey a commutation rule in order to 
comply with the symmetrization postulate: 

5(x-y) = [a x ,4]. (43) 

So far, we have introduced the quantum field a in a 
position representation a x . However, in performing ac- 
tual calculations other representations are often more fa- 
vorable, for example, bare-harmonic oscillator-, or self- 
consistent Hartree-Fock states. Consequently, we will 
delay that choice and work with a yet unspecified ba- 
sis that spans the same single-particle Hilbert space 
TL = Span({|gi) \ q\ E Q}). Here, the index set Q en- 
compasses all possible single-particle quantum numbers 
triples qi. In this generic basis, the removal of a particle 
from a position x transforms into 

Sx = 5^(x|«i)a«ft =(x|l)oi, (44) 

S qi,Q2 = [ a <m««J ■ ( 45 ) 

Here, a qi denotes a bosonic operator that removes a par- 
ticle from a general mode \qi). To simplify the notation, 
we drop the name of the dummy summation variable, 
i.e., a qi = Si, as well as the summation symbol itself. 

The temporal evolution of a weakly interacting bosonic 
gas is governed by a Hamilton operator that consists of 
a single-particle energy in the presence of the external 
trapping potential, a two-particle interaction potential, 
as well as higher order contributions. Within the dilute 
gas limit, we want to assume that these higher order con- 
tributions, mediated through three-body collisions, are 
unlikely to occur, and we disregard them. Thus, we as- 
sume a Hamilton operator of the following form: 

H = HM+HW = (0) 12 a\a 2 + «/» 1234 a\4a 3 a 4 . (46) 



Here, denotes a single-particle Hamilton operator 

with matrix elements H^ 12 = (1| p 2 /(2 m) + F ext (x) |2). 
To be specific, we assume for convenience that the ex- 
ternal trapping potential is an isotropic harmonic oscil- 
lator, i.e., K x t(x) = muj 2 (x 2 + y 2 + z 2 )/2. In most of 
the present experiments, the two-body interaction poten- 
tials Vbi n (xi — X2) are repulsive, of short range, and are 
described by the two-particle matrix elements: 

<f> 123i = <1| ® (2| tWxr - x 2 ) |3> <8> |4) , (47) 

^1234 = ^1243 = ^2134 = ^2143_ (4g) 

Only the symmetric part of the two-particle matrix el- 
ement 1234 is physically relevant. Therefore, we have 
explicitly (S) symmetrized it. In the low kinetic energy 
range that we are interested in, s-wave scattering is the 
dominant two-particle scattering event |4^ [h|. Thus, by 
discarding all details of the two-particle potential, we can 
describe the interaction strength with a single parameter 
Vq related to the scattering length a s by Vq — Ai:h 2 a s /m. 
This limit corresponds to a singular interaction potential, 
i.e., Vbin(xi,X2) = Vb<S(xi— X2). In the case of this delta 
potential, one finds for the two-body matrix elements: 

1234 = ^/ d3 x (i| x )( 2 |x)<x|3)(x|4>. (49) 
* J -00 

which need not be symmetrized, as they are symmetric 
already. However, considering the caveats that are re- 
lated to the singular functional form of the two-particle 
potential |39|| , we will only rely on the existence and sym- 
metryof the two-particle matrix elements as defined in 
Eq. @. 

It is interesting to note that the Hamilton operator 
Eq. ([hjj) is more general than its intended use. In case 
of trapped atoms that have several internal electronic 
states, it is only necessary to combine all external and 
internal quantum numbers into the definition of a single 
particle state, i.e., = |ni, l\, m\\ i*i, Mi, . . . ). This 
implies that all derived results also hold true for multi- 
component systems, provided the matrix elements -ff^i2 
and (p 1234 are generalized accordingly (for example, dou- 
ble condensate mean- field equations in Refs. p7| , p8[ ) 

The renormalization potential Qj~v accounts for the 
mean-field shifted energies that affect the single-particle 
propagation between consecutive collision events. We 
have seen already in the previous section that the mean- 
field potentials, which result from a first order calcula- 
tion Eq. (pl2|), are invariant under the particular choice 
of a energy partitioning. Now, by anticipating the re- 
sults of the first-order calculation, we do not just find a 
single mean- field potential, but several: i.e., one occur- 
ring in the equation for mean-fields, and a different one 
occurring in the equation for fluctuations. Furthermore, 
due to the inclusion of the anomalous fluctuations (see 
next section), we no longer have particle- number con- 
serving mean-field potentials either. However, as we will 



7 



use the particle-conserving part of the mean-field ener- 
gies to determine the "best" single-particle basis, we will 
also choose a renormalization potential that is number 
conserving, i.e., 



_ -1234 ~t n 23 - 



(50) 



where Q 2 ^ are the matrix-elements (yet to be deter- 
mined). 



B. The set of relevant operators 

The derivation of the kinetic equations is based on the 
premise that only a few fundamental variables determine 
the gas's evolution on a coarse-grained time scale. These 
quantities will serve as master variables and any single 
time observable is linked to them. It is an intricate ques- 
tion to determine a set of relevant variables from general 
grounds up, but we are guided by the following physical 
arguments: 

In the case of kinetic temperatures well above the tran- 
sition temperature, it is sufficient to consider only the re- 
distribution of populations f qi — (fit a qi ) within generic 
quantum levels |gi). The quantum average is defined by 
Eq. ©. 

However, as temperatures are lowered, the spatial ex- 
tension of a single-particle wave function becomes com- 
parable to the mean inter-particle distance. Thus, it will 
be necessary to consider spatial coherences as well, i.e., 
fqi >?2 = ) • Note, this is still a single-particle quan- 

tity. 

On the other hand, the most salient and fascinating 
feature of Bose-Einstein condensation is the formation 
of a macroscopic many-particle mean-field a qi = (a qi ). 
By now, this is an experimentally well established fact 
and to a large degree the mean-field is described by the 
Gross-Pitaevskii equation. Moreover, there are theoreti- 
cal predictions jl7],[49| indicating that anomalous fluctu- 
ations play a significant role as well. Consequently, we 
will also consider anomalous averages m qi ^ q2 = (a qi a q2 ) 
as independent relevant variables. 

Thus, guided by the forgoing arguments, we choose the 
set of relevant operators as 

(51) 



iel} = 


{1, 


^11 > ^92 ' 




1 'ni2 ~ 






a qi )) 


TO 9l92 = 


: {a qi 


— a qi )(a q2 


a q2 )) 


n <?l92 = 




-««i)( a » 





and denote the corresponding expectation values 7$ = 
Tr{7i cr {7} } by 



{7,|i £1} = {1, a qi , a* 2 , 

/<3l92 7 TO 9l92J n 9l92 



(52) 



With this set of independent variables, i.e., the mean 
fields and the fluctuations around them, we can parame- 
terize the reference distribution as 



a {a,cx',f,m,n} ~ eXP (°{/,m,n} 
,12 



r ^ 12 

{fi™ fa } 

l 12 * 



-mi 2 Af ? _ -ni 2 A!V_ ). (53) 

[f,m,n} {f,m,n}J v ' 

Here, we have used the implicit summation convention 
described in Eq. (ji^). The conjugate thermo-dynamic 
coordinates {O, T, A} are implicitly defined by the quan- 
tum averages 



that is 

(fa) {a,a* ,f,m,n} 



fl2, (mu){ ctja *J l m,n} = ™'W' 



(54) 



(55) 



as well as their complex conjugates n\i — m* 2 . The 
average of any other multiple operator product, occur- 
ring during the evaluation of the kinetic equations, are 
greatly simplified by the Gaufiian structure of the refer- 
ence distribution. A set of factorization rules, known as 
Wick's theorem |p7[ , can be derived and its main results 
are outlined in Appendix |b|. 



C. Renormalized master equations 

In this section we present the results of applying the ki- 
netic master equations Eq. p^ ) to the set of relevant op- 
erators defined in Eq. ( [si"| ) . Within the limits of the phys- 
ical approximations we have obtained a self-consistent set 
of equations for the mean-field amplitude a that gener- 
alizes the Gross-Pitaevskii equation, a quantum Boltz- 
mann equation for the normal fluctuations (depletion) / 
and the anomalous fluctuations rh. The large number of 
individual algebraic transformations (« 10000) that are 
necessary to obtain the final result prohibits attempts 
to evaluate the collision terms manually. Therefore, we 
developed a symbolic algebra package that performs the 
required calculations. The presentation of the final re- 
sults of this calculation is greatly simplified by introduc- 
ing the following single-particle Hilbert-space vectors (co- 
, contra- variant) 



<7i,<?2 G Q}- 



(a) = \a) =a x |1) , (fi)t = ( a \ = a\ (1| 
normal operators [tensor rank (1,1)] 

/=/ 12 |1)(2|, /W=a5ai |1)(2|, 
pseudo operators [tensor rank (2,0)] 

m=m 12 |1)|2), m^=a 2ai |1) |2) , 

and their Hcrmitian conjugates n = ff, 



(56) 
(57) 
(58) 
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1. Mean fields 

For the macroscopic many-particle mean field \ot), we 
find 

j t \a) =-i(H^ + lU f ^+2U f -)\a)+ (59) 
-iVfnZ. {a\+Lf ] . - [a], 

with a collision term given by 

£ {L,a-,/,m,n}® = ( T ffU+f) ~ T (^+f) (1+/) / + ( 60 ) 

+ 2r /rraS -2r (1+/)fflS ) |a) + 

+ 2 ( r /m(l + /) - r (l + /)m/) Z H ■ 

First of all, the field evolves unitarily in the nonlinear 
hermitian Hamilton operator which consists of a free part 
H(°\ determined by the external trapping potential and 
any other applied electro-magnetic fields. The second 
and third contribution are the collision induced mean 
field potentials, denoted by Uf<. C ) and Uj. While the first 
of these potentials is proportional the mean field den- 
sity itself, the second potential Uj reflects the influence 
of the normal fluctuation upon the mean field. It is im- 
portant to note the different weighting factors 1 and 2 
multiplying the potentials. They arise from the different 
quantum statistical fluctuation properties of a c-number 
mean field and a normal single particle density. Exactly 
the same weighting factors are also found with the varia- 
tional Hartree-Fock-Bogoliubov approach [|l7j • The mean 
field potential is defined in terms of the two-particle in- 
teraction matrix elements Eq. ( [47| ) and a single-particle 
density operator / that can be either or / 

C/ / = 2 12 ' 3 ' 4 '/ 3 '2' |1)(4'|. (61) 

Due to the Hermiticity of the two-particle interaction en- 
ergy and the positivity of the single particle density, it is 
also self-adjoint [j], i.e., Uf — Uj. It is interesting to see 
that in the case of a delta-potential Eq. ([l9|) and a scalar 
field (i.e., no internal degrees of freedom), the mean field 
potential reduces to the well known potential energy den- 
sity that is proportional to the local mean field density 

(x|[/ /(c) |y) =V^(x-y)|(x|a)| 2 . (62) 

The fourth linear collisional contribution is proportional 
to the anomalous coupling strength Vfh contracted (Z) 
with the adjoint field (a\. In general, we obtain the con- 
traction of two tensor fields A/LB, from a basis represen- 
tation of the two fields and a subsequent contraction of 
the last index of A with the last index of B. 



^This potential Uf is not related to the single-particle prop- 
agator fjfl (r) defined in Eq. @ 



This non-Hermitian coupling is in general mediated by 
an anomalous average m, and explicitly given by 

V m =2c/> 12 ' 3 ' 4 ' m 3 , 4 , |1)|2'). (63) 

In here, m stands for any anomalous average, either m^°' 
or rh. From the definition of the anomalous coupling, it 
can be seen easily that V m = is symmetric. 

All of the remaining terms in Eq. (|6^) are second order 
collisional contributions. They always appear in pairs 
where one term corresponds to an in-process while the 
sign reversed companion describes a loss out of the field. 
A closer inspection reveals that there are essentially four 
types of processes occurring, i.e., collision events that in- 
volve zero to three anomalous averages. Furthermore, 
one finds that a normal fluctuation / on the in-side is al- 
ways accompanied by a bosonically enhanced (1 + /) on 
the out-side, and vice versa. On the other hand, when- 
ever a mean field density f( c \ an anomalous mean field 
density m^ c \ or an anomalous fluctuation rh occurs in 
an in-process they appear unaltered on the out-process. 
This behavior is analogous to atomic transition rates de- 
scribed by the Einstein A- and B coefficients which can 
be attributed to stimulated absorption- and emission, as 
well as spontaneous emission processes. The fact that the 
mean-field is never bosonically enhanced supports the in- 
terpretation that the mean-field acts as a classical driving 
field. 

In detail, these collisions operators are described by 
the following operators and pseudo-operators 

Tfff = 8^ 2 ' 3 ' 4 '<' 2 " 3 " 4 "/3'l"/4'2"/4"2' |1> (3"| , 

T fmf = 8^ 2 ' 3 ' 4 >;" 2 " 3 " 4 7 3 'i»m 4 -3«/4»2' |1> |2") , 
T fmn =8<b 12 ' 3 ' 4 'rt' 2 " 3 " 4 "f3'vm i , 3 „n 2 „ 2 , |1) (4"| , 
T mmn =8^ 2 ' 3 ' 4 >;" 2 " 3 " 4 "m 3 - 4 ,m wtl2 » 2 - |1) . 

(64) 

From the time average over the interaction picture 
Hamilton operator that appears in the kinetic equation 
Eq. (|4^), one obtains an approximately energy conserv- 
ing two-particle matrix element. 

^' 2 " 3 " i "(t) = f dre"V 234 x (65) 

J — OO 

^ { t 7W}( T ) 1 "x iir {7W}( T ) 2 " 2 ^7(t)}W 3 '>{7(t)}^) 4 "4' 

The restricted propagator that has been used here is ex- 
plicitly given by the time-ordered exponential 

K Mt)}{r) =Te'J'°' is ( H(D)+g ft^W'»»). (66) 

where we obtained the operators and Q{ 7 } from the 
matrix elements given in Eqs. ((46|,j50|) . 

To see qualitatively why <f>* 2 3 4 is essentially non- 
zero only on the energy-shell of thickness rj, it is useful 
to represent the restricted propagator with respect to the 
eigen-states of 
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(H^+Q {im )\l) = £l (t)|l). (67) 

By assuming that the energy levels change adiabatically 
slow, one obtains approximately 

j.l"2"3"4" _ /l"2"3"4" ' - ' ■ - — ' 



'2"3"4" 

(68) 



which is non-zero only if the energy difference 

Ai" 2 "3"4" = + £2"{t) - £3"(t) - £4"(*) is smaller 

than r\. 

n %^J K =-^(A) + iV^, (69) 

This result is analogous to the second order Born-Markov 
approximation 6a|. 
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2. Normal fluctuations 



The kinetic equation for the normal fluctuations (de- 
pletion) generalizes the quantum-Boltzmann equation 
found in many textbooks PJ37|] 

jj = ~i [H {0) + 2 U fM + 2 Uj, f] + (70) 

—i V( m (,c) . ~,i Zn+imZy/ ,„,,-,+£, , ~ c J/1, 

l m '-Tin) (■m<. c >+rn) {a ,cc* ,f ,m ,n} lJ J ' 

with second order collisional contributions 

L {a,a*,/,m,n} ^ = ( r //(l+/) + 2 /(*+/) + r ///W + 2 T f(.mM+m)n + 2 + 2 r /Wmn) (* + /)+ ( 71 ) 

~( r (l + /)(l + /)/ + 2r / (c) (l+/)/ + r (l+/)(l + /)/ (c) + 2r (l+/)(m<'=)+m)S + 2r (l+/)m»W + 2r /Wfts)/ + 
+ 2 ( r /(mW+m)(l+/) + r /<'=)A(l+/) + r /m/<<=) _ r (l+/)(wW+m)/ ~ r /(<=)m/ _ r (l + /)m/(<=>) ^ « + h.C. 

First, one finds a unitary evolution in the presence of 
the external trapping-, the mean-field-, and the normal 
potential. Both of the self-induced potentials, f7y(c) and 
Us are weighted by a common factor of 2 (compare ) . 
This is in contrast to the weighting factors appearing the 
mean- field Hamilton operator Eq. (f>9|). But again, this 
fact can be traced back to different quantum statistical 
properties of the mean-field and the fluctuations. Sec- 
ond, it can be seen that the anomalous coupling strength 
is now proportional to the total anomalous average, i.e., 
mA c ) + fh. Third, in the absence of any mean-fields 
or anomalous averages the second order contribution in 
Eq ( |7l| ) reduces to the well known Boltzmann collision 
term 

r /7(i+/)( 1 + f) ~ r (i+/)(i+/)/^ ( 72 ) 

By further assuming that the normal fluctuations are pre- 
dominantly diagonal in an energy eigen-basis defined by 
the non-linear Hamilton operator of Eq. ( |70| ) (ergodic hy- 
pothesis), one recovers the Bose-Einstein distribution as 
the stationary distribution of particles within the quan- 
tum levels. However, the presence of the mean- field, as 
well as the anomalous averages lead to additional collision 
processes that must not be ignored in general. Eventu- 
ally, these processes will lead to a self-consistent equilib- 
rium partition of particles between mean-fields, normal- 
and anomalous fluctuations. While a detailed numerical 
self-consistent solution of the set of kinetic equations is 
still under investigation, it is important to see that the 
total particle number (N) = Tr{/( c )} + Tr{/} is always 
conserved [compare Eq. (J40|)] 

|(iV)=0. (73) 
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3. Anomalous fluctuations 



In contrast to the normal fluctuations, the anomalous 
fluctuations do not evolve unitarily but rather as a tensor 
of rank (2,0). Both, left and right generators of the time- 
evolution are identical to the Hamilton operator of the 
normal fluctuations. 

-|m = -i (H {0) + 2 U fM +2U f ) Zm-imZ (H {0) + 2 C/ /(c) + 2 U f ) + (74) 
-*V( m (<o +ffi) Z(l + /) -ifZV {m(c)+fh) + L fl a ,j^ } [m}, 

£ S,o',/,m,S}[ ft ] = ( r /7(l+/) + r ///«=> +2r //W(l+/) + 2r /m(n«=)+n) + 2T fm(On) Z ™+ ( 75 ) 
_ ( r (l+/)(l+/)/ + r (l+/)(l+/)/< c ' + 2r (l+/)/(<=)/ + 2r (l+/)m(nW+n) + 2 T (l+f) m Mn) Zm + 
+ ( 2r /(m( c )+ft)(l+/) + 2r /m/(=) + 2r /( c )m(l+/) + r mro(n(<=)+n) + 2 r ftm (c )? ^ Zf + 

-( 2r (i+/)(mW+m)/ + 2r (i+/>/W + 2r /(c)™/ + r A f fl( „( c)+fl) +2r ftmWs ) ^(1 + /) +transp. 

These three sets of master equations for the mean- 
field, the normal-, anomalous fluctuations constitute the 
main result of this article. They unify and generalize 
simpler equations which have been obtained previously 
also by other methods. However, so far, we have not 
discussed the physical implications that will arise from 
a self-consistent solution of these equations. Specifically, 
we need to determine the following problems: (I) the 
equilibrium distribution of particles partitioned between 
mean- field, normal-, and anomalous fluctuations; (II) the 
importance and quantitative size of anomalous fluctua- 
tions; (III) the collisional damping rates and second order 
energy-shifts; (IV) the response of the equilibrium system 
to weak external perturbations, i.e. the collective exci- 
tation frequencies via linear response theory; (V) critical 
phenomena occurring around the onset of condensation; 
or, for example, (VI) the dynamics of the growth of the 
condensed phase. 
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IV. OUTLOOK 

In the previous section, we have enumerated several 
quantities that need to be determined and interesting 
paths along which detailed calculations could proceed. 
We believe that amongst these issues, it will be most 
crucial to address problems (I) and (II) around T ~ 0, 
first. On one hand, present-day experiments have estab- 
lished that the mean-field description yields good agree- 
ment. On the other hand, there are various approaches 
to the self-consistent equilibrium for normal and anoma- 
lous fluctuations and not all implications have been elu- 
cidated. 

The standard route to investigate this problem is based 
on finite temperature calculations in the Hartree-Fock- 
Bogoliubov description. Various schemes employing, for 
example, the quasi-static Popov approximation, or more 
dynamical methods that go beyond it (i.e. the collision- 
less regime) are being investigated by several research 
groups. 

This present, non-equilibrium approach provides an al- 
ternate route to the stationary solution. In particular, we 
expect that the presence of a large condensed phase will 
lead to a strong correlation of the low energy part of the 
normal and anomalous fluctuations (ss 1 — 2 times the 
chemical potential \i of the condensate), while the high 
energy tail will be mostly in detailed balance at some 
temperature T. However, such a macroscopic "polariza- 
tion" of the low energy part of the fluctuations can not be 
described within a simple ergodic hypothesis, and there- 
fore requires a full quantum treatment. 

The main obstacle to overcome in numerically answer- 
ing this problem is the unfavorable scaling law of the 
collision operators. From a simple operations count, one 
finds that there are TV 8 summations involved if N is the 
number of energy levels being considered. 

This burden can be alleviated by being more specific, 
i.e. by postulating a completely isotropic situation for 
a single condensed phase, an isotropic trapping poten- 
tial, a rotationally invariant initial condition, as well as a 
short-range central two-particle interaction. Within this 
simplified model, one can then decompose all involved op- 
erators in terms of angular momentum sub-manifolds (i.e. 
irreducible tensor sets and use of Wiger-Eckart theorem) . 
This assumption makes the quantum mechanical treat- 
ment of the low energy region (fi < e < 2 /i, N « 10 — 20) 
feasible and will lead to a self-consistent equilibrium. A 
detailed numerical investigation is in progress and results 
will be reported. 



CONCLUSIONS 

In this article, we have revisited the Chapman-Enskog- 
Bogoliubov procedure of non-equilibrium statistical me- 
chanics to describe the kinetic evolution of a condensed 
bosonic gas of atoms towards equilibrium. Within a 



second order Born-Markov approximation, we consider 
the collision dynamics of macroscopic mean-fields, nor- 
mal fluctuations, and anomalous averages. In particu- 
lar, we have obtained a coupled set of master equations 
for these quantities that encompass the Gross-Pitaevskii 
mean-field equation, as well as the quantum-Boltzmann 
equation for the normal fluctuations as limiting cases. 
The mean-field potentials that are obtained from a first 
order calculation are in agreement with the results of 
a variational Hartree-Fock-Bogoliubov calculation. Be- 
yond these first order energy shifts, we obtain second 
order collisional energy shifts and damping rates that 
are bosonically enhanced. We expect our results to be 
valid when strong collisions are well separated in time 
and when the mean-field induced energy shifts m ay b e 
neglected during a strong collision event [see Eqs. JA7|), 
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APPENDIX A: 

With the definitions for the interaction-picture rep- 
resentation, we can rewrite the commutator term of 
Eq. (||) as 

{ 7 }V / {7(t;{7})}' {7(t;{7})} j ' 



(Al) 



The trace term evaluated along the trajectory simplifies 
to 



Tr 



{^S(TS{7})}' Ti] Cr {^(r;{7})}} _ ^ 



Finally, the self-tuning term of the reference distribution 
gives 



(A3) 



fry, (t; {7} A 



where we introduce auxiliary operator valued vectors 
D{ 7 }(t) and matrix- valued coefficients, Sm(t) by 



(0)V 



<S {7 }(t) = I — J K {i} {t)i . 



(A4) 
(A5) 
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If these results are put together, one finds for the first 
order correction of the coarse-grained statistical operator 



'{7} 



dr 



+ {[D M (r) u a™] +d^\) 5 {7} (r)« 
Tr{[H^ } (r),^]a^ } }). (A6) 

This expression is formally equivalent to Eq. (|2^) . How- 
ever, a closer inspection of the S and D terms shows that 
they contain higher order energy corrections induced by 
renormalization energy Q/ 7 \- In particular, from a short 
time Taylor-expansion one finds that 



5 {7} (t) = 0-iO [r d lz Q {1) 



{7} 



(A7) 
(A8) 



Consequently, we will disregard the effect of the mean- 
field onto the temporal evolution of -D{ 7 }(r) and S*{ 7 }(r) 
during a strong collision event and replace them by their 
"bare" values attained in the absence of the mean-field 
shift. 



APPENDIX B: A GENERALIZED WICK'S 
THEOREM 

The Gaufiian structure of the reference distribution 
Eq. (|5^) is particularly useful, as it permits the system- 
atic application of Wick's theorem p7j . This is a set of 
rules to efficiently evaluate quantum averages for multi- 
ple operator products as 



(Bl) 



In this average, for example, the operator ipi represents 
either an operator a± or a\. 

First, the displacement rule shifts any operator ipi by 
its c-number expectation value "01 which is either a,\ or 
a\ , and replaces the quantum average by an average that 
has zero mean values: 

= ((V>1 + i>l)(i>2 + ^2) • ■ ■ {i>n + V'0){o,0,/,fn,n}- 

Second, after expanding the multiple products, one can 
discard all averages that involve an odd numbers of op- 
erators: 



(B3) 



And third, for the remaining averages, one can use the 
Gaufiian factorization rule: 



<lM2---^>{0,0,/,f&,n} = (B4) 
= (i>1^2}{0,0,f,ih,n}(^3 ■ ■■^){0fi,f,fn,n} + 
+ (^1^3){0,0 ) /,ft,a}^^*'--^ 2 '){0,0 1 /,ln,n} + 

+ (&$2.){o,0,/,iii,n}v& ■ ••^2.-l){0,0,/,a,5i}- 

By proceeding recursively, one has finall y ev aluated the 
complete multiple operator average Eq. (Bl). 
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